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ABSTRACT 

The cold dark matter theory predicts triaxial dark matter haloes. The radial distribu- 
tion of halo ellipticity depends on baryonic processes and the nature of dark matter 
particles (coUisionless or collisional). Here we show that we can use lensing flexion 
ratios to measure the halo ellipticity as a function of radius. We introduce a weight 
function and study the relationship between the first and second order statistics of 
flexion ratios, both of which can be used to reduce the bias in the estimate of ellipticity. 
We perform numerical tests for our method, and demonstrate that it can reduce the 
bias and determine the halo ellipticity as a function of radius. We also point out that 
the minimum mean flexion ratio can be used to trace the centres of galaxy clusters. 
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1 INTRODUCTION 

Structures of cluster- and galactic-sized dark matter haloes 
predicted by N-body simulations of Cold Dark Matter 
(CDM) have several important features, e.g. a nearly uni- 
versal radial density profile for haloes over a wide range 
of masses (Navarro et al. 1996, 1997), and triaxial shapes 
(Jing & Suto 2002; Springel et al. 2004). These features are 
related to the nature of dark matter particles, the forma- 
tion processes of galaxies and clusters, and their environ- 
ments (Kuhlen et al. 2007; Maulbetsch et al. 2007; Debat- 
tista et al. 2008; Faltenbacher et al. 2009; Bett et al. 2010; 
Ragone-Figueroa et al. 2010; Wang et al. 2010; Vera-Ciro 
et al. 2011). For example, major mergers of dark matter 
haloes may play an important role in forming their shapes 
(van Haarlem & van de Weygaert 1993). Mismatches be- 
tween luminous galaxy shapes and total mass shapes provide 
evidence for the existence of dark matter (Suyu & Halkola 
2010). Numerical simulations under different assumptions 
predict dark matter haloes with different shapes (Bullock 
2002; Bailin & Steinmetz 2005). For instance, coUisionless 
CDM predicts prolate triaxial haloes (Allgood et al. 2006), 
while simulations with gas cooling suggest more spherical 
shapes of dark matter haloes (Kazantzidis et al. 2004) . More- 
over, taking into account the baryonic cooling effect, inner 
halo shapes are often found to be more spherical than their 
outer counterparts (Springel et al. 2004; Kazantzidis et al. 
2004). In simulations without baryonic cooling, however, the 
opposite trend has been seen, i.e. the outer parts of haloes 
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appear to be rounder than inner parts (Dubinski & Carlberg 
1991). 

The shapes of cluster-sized haloes or galaxy-sized haloes 
can be studied with different approaches, such as X-ray 
observations and the Sunyaev-Zeldovich effect (Reblinsky 
2000; Lee & Suto 2004; Wang & Fan 2004; De Filippis et al. 
2005; Sereno et al. 2006). One can also use the spatial dis- 
tribution of satellite galaxies as a tracer for the shapes of 
their host dark matter haloes (Wang et al. 2008, 2010). 

Gravitational lensing provides a powerful tool to study 
the mass distributions of clusters of galaxies as well as galaxy 
haloes (see Bartelmann & Schneider 2001; Refregier 2003; 
Schneider et al. 2006; Munshi et al. 2008, for reviews on 
weak lensing). This is because gravitational lensing probes 
the distribution of matter regardless of whether it is lumi- 
nous or dark. Weak lensing techniques have been used for 
cluster mass reconstructions (e.g. Clowe et al. 2006; Bradac 
et al. 2008) and for determining halo ellipticities (Hoekstra 
et al. 2004; Corless et al. 2009; Deb et al. 2010; Howell & 
Brainerd 2010). The method of stacking galaxies to deter- 
mine galaxy halo ellipticity was proposed in Natarajan & 
Refregier (2000) and has been used to determine ellipticities 
of both cluster- and galactic-sized haloes (Evans & Bridle 
2009; Mandelbaum et al. 2006). For example, a mean ellip- 
ticity of 0.46 has been found from a sample of 25 clusters 
using Subaru data (Oguri et al. 2010). 

Flexion has been studied as the derivative of the sur- 
face mass density and the shear, responds to small-scale 
variations in the projected mass distribution (Goldberg & 
Natarajan 2002; Goldberg & Bacon 2005; Bacon et al. 2006). 
Different techniques have been developed to measure flexion 
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(Irwin & Shmakova 2006; Okura et al. 2007; Massey et al. 
2007; Schneider & Er 2008; Velander et al. 2010; Fluke & 
Lasky 2011; Cain et al. 2011). Due to the complexity of 
deconvolving for the point spread function, current lensing 
image measurements are not accurate enough to fully study 
flexion. It has only been implemented on limited samples 
of real data. Okura et al. (2008) and Leonard et al. (2007) 
have performed flexion measurements to detect substruc- 
tures in the cluster Abell 1689 using Subaru and HST data. 
Though flexion noise is still poorly understood, flexion can 
potentially contribute to cosmology in several aspects, such 
as exploring dark matter haloes of galaxies and clusters, 
especially substructures (Leonard et al. 2009; Bacon et al. 
2010; Er et al. 2010), understanding the large scale struc- 
tures through cosmic flexion (Munshi et al. 2011; Schaefer 
et al. 2011). Hilbert et al. (2011) also proposed to reduce 
the distance measurement errors of standard candles using 
lensing shear and flexion maps. In Hawken & Bridle (2009), 
galaxy halo ellipticity has been studied with flexion. Re- 
cently, Er & Schneider (2011, hereafter ESll) developed a 
new approach using the ratio of the tangential to radial flex- 
ion to estimate the halo ellipticity as a whole. The approach 
is independent of lens strength, i.e. the mass or the lens 
redshift. In particular, the flexion ratio is independent of 
the halo-centric distance. In this paper, for the first time, 
we use the flexion ratio to estimate how the halo ellipticity 
varies with radius. 

In reality, there are several complications. For example, 
the background galaxy number density may not be suffi- 
ciently high, and galaxies may not be circularly symmetri- 
cally distributed. Even worse the intrinsic and other noises 
will cause an overestimate for the mass ellipticity. In this pa- 
per, we introduce a weight function and a method that com- 
bines the first and second order statistics of flexion ratios to 
reduce the bias (see Appendix A), we test our method with 
numerical simulations (Sect. 3). Our results are discussed in 
Sect. 4. The cosmology that we adopt throughout this paper 
is a ACDM model with Q.a = 0.75, Vim = 0.25, ag = 0.9 and 
a Hubble constant Hq — 73kms^^Mpc^^. 



2 BASIC FORMALISM 

The full formalism described here can be found in Bacon 
et al. (2006); Schneider & Er (2008). Weak lensing shear 
and flexion are conveniently described using a complex for- 
malism. We adopt the thin lens approximation, assuming 
that the lensing mass distribution is projected onto a sin- 
gle lensing plane. The dimensionless projected mass density 
(convergence) can be written as k,(0) — E(0)/Ecr, where is 
the position coordinate, E(0) is the projected mass density 
and Ecr is the critical density, given by 
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for a fiducial source located at redshift Zs . Here Ds , Dd and 
Dda are the angular diameter distance between the observer 
and the source, the observer and the lens and between the 
lens and the source, respectively. 

The first-order image distortion is measured by shear 7, 
which transforms a circular source into an elliptical one. The 
second order effect, which is called fiexion, is described by 



two parameters. The spin-1 fiexion is the complex derivative 
of K 
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and the spin-3 fiexion is the complex derivative of 7 

g = Vc7- (3) 

2.1 Flexion ratio 

The full definition of the flexion ratio can be found in ESll, 
here we only briefly describe the essential quantities. The 
flexion ratio is given by 
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where J^t and J^r are the tangential and radial components 
of the spin-1 flexion J^, respectively. ESll has shown that 
for an elliptical isodensity contour of the mass distribution, 
the flexion ratio provides an elegant estimator for the halo 
mass ellipticity e and the orientation angle (j>o 
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where is the polar angle, e = {a — b)/{a + b) and a, b stands 
for the major and minor axes of the projected halo ellipse 
respectively. The expectation value of the flexion ratio is 
thus directly related to the halo ellipticity 
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We then have two estimators e and e for the halo ellipticity 
using (r) and {r^) 
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Due to the intrinsic and other noises, both estimators are 
biased. The bias behaviors in (r) and (r^) are different. By 
comparing (r) and (r'^), we can obtain an approximation 
for the bias br (Eq.All), which we use to subtract from the 
observed (r^ts) by 



(»"dcan) = (^obs) - 2{r)br 



(11) 



Notice the bias br is not zero since it is not the same as 
the average noise, which is zero. Using (rciean) in Eq.(lO), a 
"clean" estimator is obtained, and will be employed in the 
following section. A detailed description of bias reduction is 
given in Appendix A2. 
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3 MASS ELLIPTICITY AS A FUNCTION OF 
RADIUS 

In this section we present how we estimate mass ellipticities 
at different radii using numerical simulations. For an ellip- 
tical isodensity contour mass distribution, the flexion ratio 
r is independent of the angular distance to the centre of 
the lens 6. Therefore the ellipticity as a function of 9 can 
be estimated. Several annular bins could be easily applied 
to the data. However, at small radii, the number of back- 
ground galaxies in each bin is low, since it is proportional to 
9 for a fixed bin width. On the other hand, the flexion signal 
drops rapidly at large 9, oc 9~^ in an isothermal halo, and 



even faster in an NFW halo, approximately oc 9 



>rs, 



where r^ is the scaling radius). One might expect a better 
signal-to-noise (S/N) ratio at smaller radii. In reality, both 
are biased for different reasons. At small radii we obtain an 
underestimate due to the lack of background galaxies, and at 
large radii the estimate is strongly affected by noise, which 
results in overestimates. 

Shapes of dark matter haloes have been studied by 
many authors using numerical simulations (Jing & Suto 
2002; Bullock 2002; Springel et al. 2004; AUgood et al. 2006). 
In these studies, there are different ways of modeling haloes 
as ellipsoids. One popular method uses some form of eigen- 
vectors of the inertia tensor. Another method fits the iso- 
density surfaces of the halo as ellipsoids (Jing & Suto 2002). 
In this paper, we employ the latter method. We fit the nth 
isodensity contour using the nth surface mass density k„ 
and estimate the major axis an and minor axis bn- We then 
calculate the ellipticity using (a„ — hn)/[an + bn) at the 
"effective" radius \/anbn- The ellipticity derived using the 
flexion ratio approach will be compared with the ellipticity 
calculated in this way. 



3.1 A halo modelled by an NFW profile 

We first generate an a halo modelled by the NFW (Navarro 
et al. 1996, 1997) profile with constant ellipticity. The di- 
mensionless surface mass density is written as (Bartelmann 
1996; Bacon et al. 2006) 
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where x is the radius in units of the scaling radius Ts, such 
that X = 9/r, and 9 = ^9l/{l + ef + 9l/{l ~ eY . The el- 
lipticity we choose is e = 0.15. We define Ks = pcrit»'sAc/Ecr, 
where Ac is the dimensionless characteristic density (see Ap- 
pendix in Navarro et al. 1997). The lens halo is at red- 
shift Zd = 0.6 and sources at Zs — 1.48. The halo has 
mass M200 = 1.8 X lO^*/i~^M0 and concentration parameter 
c — 7.2 (the corresponding scaling radius rs — 27 arcsec). 
The function f{x) is given by 
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The fiexion data are generated by 
2J-, 
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where J^s ~ KsDd/rs- The polar angle can be determined 
from tan (/> = -rz ,„ „ , and 
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Two sets of mock data are generated. First we generate flex- 
ion data on 100 x 100 grids in a field of 2 x 2 arcmin^. In 
the other data set we have 1000 flexion data on random 
locations in the same field. We add the intrinsic Gaussian 
fiexion noises 0.03 arcsec"^ for the 9i and 9^ components 
respectively (more details are given in Appendix B). The 
data are discarded if they are located too close to the halo 
centre, 8 < Omin ~ 6 arcsec or too far from the halo centre, 
9 > 6'max ~ 55 arcsec. Since a very large flexion cannot be 
measured (Schneider & Er 2008), data points with ] J^| > 0.5 
are also discarded. Moreover, due to the upper bound of flex- 
ion ratio (ESll), those with r > 4.5 are discarded as well 
(in all the numerical tests below, similar fllters but different 
Smax, duiin are employed). The ellipticity e is also estimated 
using the surface mass density k as a comparison. We apply 
8 bins in performing the flexion ratio method. In each annu- 
lus, an estimate for e is obtained, as well as the radius, taken 
as the mean radius over all the data in the annulus. The re- 
sult is shown in Fig. 1. For data with noise, 50 realisations 
are used. The mean values and standard deviations of our es- 
timates over 50 realisations are shown by the diamonds and 
error bars. One can see that the ellipticity estimates using k 
(filled circles) and grid noise-free data (crosses) closely trace 
our input values. The results using data with noise have a 
mean error < 10% at small radii (< 30 arcsec). However for 
the large radii, the results are dominated by the intrinsic 
noise, and thus are significantly overestimated as expected. 



3.2 A galactic-sized iialo 

We take a galaxy halo from the Aquarius project (Springel 
et al. 2008), which is a suite of A'^-body simulations of galaxy- 
sized CDM haloes. The halo Aq-F at resolution level 2 is se- 
lected for our test, with a softening length of ~ 0.05h~^ kpc, 
mass resolution of 5 x 10''/i~^Mq. The main halo contains 
about 1.4 X 10* particles - a total mass of 7.2 x 10^^ ^"^M© 
within its virial radius of 153/i~^ kpc. We place the halo at 
a lens redshift Zd = 0.6. The halo is then projected along 
20 different directions to generate 20 different sets of mock 
data. The background sources are randomly located in a 
"source plane" at redshift Zs = 1.0. The k and flexion maps 
are generated on 100 x 100 grids in the lens plane. The flex- 
ion for each image is linearly interpolated using values of 
the four nearest grid points. Furthermore, we add Gaussian 
intrinsic flexion noise (0.03 arcsec"'^) to the data. 

To test the ellipticity as a function of radius, we perform 
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Figure 1. Comparison of different halo ellipticity estimators us- 
ing an NFW halo. The filled circles are estimated directly from 
fitting with surface mass eUipse, the crosses are estimated from 
the flexion ratio method (Sec. 3.2) with noise- free grid-based data. 
The diamonds arc estimated using the clean method but with 
noisy data, and they are mean result over 50 realisations, the 
errors represent the standard deviations. 



our method on stacked galaxy-galaxy tensing since tlie num- 
ber density of background galaxies is not sufficient to study 
an individual galaxy. In each mock data projection, we use 
10 flexion at random locations in an area of 20 x 20 arcsec^ 
(which is equal to 90 background galaxies/arcmin^). The 
data from 20 different projections are aligned and stacked 
with repect to the centre of the lens halo. 50 realisations are 
generated to calculate the mean value. Similar filters as in 
the last section are employed here, except that we discard 
data closer than 1 arcsec o the halo centre, or further than 8 
arcsec away. The ellipticity using the surface mass density is 
calculated, and the average ellipticity over 20 projections is 
shown by tire filled circles in Fig. 2. The open (blue) circles 
are the estimates from the flexion ratio method using data 
without intrinsic noise and the diamonds (purple) are those 
using noisy data. They are both averages over 50 realisa- 
tions. The error bars represent tlie standard deviations in 
each bin. The error bars using noise free data in Fig. 2 are 
caused by the variances of halo ellipticities projected from 
different directions, while the error bars using noisy data 
are affected by the estimate error. We can see that the error 
bars increase with radius, as intrinsic noise becomes more 
important. 



3.3 A cluster-sized halo 

In addition, we test our approach on an A''-body simu- 
lated galaxy cluster halo. For galaxy clusters, we liave more 
background galaxies in "observations". Tlierefore the flex- 
ion number density is high enough to resolve an individ- 
ual cluster. The test cluster halo PII-A-2 is taken from 
the PHOENIX simulations (Gao et al. 2011 in prep). At 
the second resolution level, the simulation has a softening 
length of 0.32/i"^ kpc, mass resolution 5 x IQp MqH'^ . The 
main halo contains 1.5 x 10* particles, with a total mass of 
6.6 x IO^Mq/i"^ within a virial radius of lAh~^ Mpc. Two 
simulated cluster haloes are used to generate mock data: one 
has only the main halo, while the other one in addition in- 
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Figure 2. Similar to Fig. 1, but for results from simulated galaxy- 
galaxy lensing. The circles are estimated using the flexion ratio 
method (Eq.6) using noise free data, and the diamonds are esti- 
mated using noisy data. The errors represent the standard devi- 
ations over the 50 realisations. 



eludes massive subhaloes that are likely to have noticeable 
signatures. The number of subhaloes is about 5% of the total 
subhaloes with mass M > 10^ Mq that are identified from 
the simulations. The halo is located at Zd — 0.5, while the 
background sources are chosen to be at Zs = 1.0. We project 
the same halo along 3 different directions to generate mock 
data. 

In Fig. 3, the surface mass density k map from one pro- 
jection (for the right panel in Fig. 4) is shown. One can see 
that there are some small scale fluctuations due to substruc- 
ture. Moreover, the shape of the halo is not perfectly ellip- 
tical. It is easy to see that there is a bump on the left side 
of the halo. Both of these will cause some difficulties in esti- 
mating the halo ellipticity since fiexion is sensitive to small 
scale variations. In each direction, we have 2 sets of data 
generated from the cluster haloes with and without sub- 
haloes. We have 50 realisations with each having 800 back- 
ground fiexion data points on an area of 200 x 200 arcsec'^ . 
Such a configuration is available via the current surveys, 
e.g. the COSMOS survey has a background density of 76 
galaxies/arcmin^ (Schrabback et al. 2010). Similar filters as 
in previous sections with 6'min ~ 10 arcsec and 6'max = 70 
arcsec are applied here. About 300 data points are obtained 
on average in each realisation. We make 6 bins in each real- 
isation, and estimate the mean ellipticity and its standard 
deviation in each bin over the 50 realisations. The results 
are shown in Fig. 4. The filled circles are directly estimated 
from the surface mass density without substructure. The di- 
amonds are the results from the data with the main halo, 
and the crosses are the results from the data with the main 
halo and the substructures. We can see that our estimates 
are slightly biased. The large error in the middle panel is 
caused by the intrinsic noise, which has very small elliptic- 
ity. In an additional test, we perform our estimates with a 
background galaxy number density of 45 galaxies/arcmin'^ 
and obtain a similar result but with slightly larger errorbars 
(~ 10%). 

There are however some other reasons that can cause 
bias in our estimate. First of all, we need to point out that 
the estimate directly from mass is not the "true value" of 
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the ellipticity. When we fit the iso-density contours using k, 
we notice that the orientations of different annuli are not 
aligned. Moreover, the annuli are not concentric, especially 
for clusters. Although our estimate using the flexion ratio 
does not depend on the radius or orientations, the binning 
of the data will be affected by the rapid variation of el- 
lipticity and orientation. Secondly, the intrinsic noise and 
substructures will usually cause an overestimate (ESll). As 
pointed out by Bacon et al. (2010), the galaxy-galaxy lens- 
ing flexion varies due to substructures. As an additional test, 
we calculate the standard deviations of tangential flexion in 
each bin for our NFW halo and the PHOENIX halo with 
substructures. In Fig. 5, the values indicated by pluses are 
significantly smaller than the others since the NFW halo 
has regular shape and contains no substructures. The oth- 
ers from the PHOENIX halo have higher deviation as we 
expect. We can also see that the ellipticity also affects the 
deviations. 

The centres of galaxy clusters are often not well deter- 
mined from observations, e.g. it is often assumed that the 
mass centroid coincides with the brightest cluster galaxy, 
which may not be correct in all cases. Also mass and light 
peaks are often offset in the weak lensing shear method (Di- 
etrich et al. 2011). On the other hand, in ESll it has been 
found that a centroid offset to the true mass centre of the 
halo causes a slight overestimate of the halo ellipticity. The 
reason is that the flexion ratio is defined with respect to the 
true centre of halo mass, a centroid offset will introduce ex- 
tra asymmetry in the estimate. Therefore, one can use the 
flexion ratio as a probe for the cluster centre, since (r) , {r^^^^) 
or an estimate of ellipticity, in principle, will reach its min- 
imum when the flexion ratio is calculated with respect to 
the true mass centre. However the intrinsic and other kind 
of noises can cause some complications. In order to see this, 
we perform another test for the centroid offset problem. We 
generate 10 realisations of mock data. For each realisation, 
we use 100 random mock centres within 20 arcsec to the 
true mass centre. The flexion ratio and (rots) ^re calculated 
with respect to the mock centre, and the result is shown 
in Fig. 6. One can see that although (r^^s) does not reach 
minimum at the true centre of the cluster, the variance of 
(^obs) monotonically increases with the offset. This is also 
because, as we mentioned in the last paragraph, the mass 
distribution of galaxy clusters is complicated, e.g. the mass 
centre may vary with the radius. Nevertheless, the minimum 
flexion ratio statistics provides a complementary approach 
for cluster centre determination. 



4 CONCLUSIONS AND DISCUSSIONS 

In this paper, we have studied the flexion ratio method for 
estimating shapes of elliptical dark matter haloes. We per- 
formed numerical simulations to study our ellipticity esti- 
mators and obtained the halo ellipticity as a function of 
radius. In particular, for galactic-sized haloes, we are able 
to measure the mean variation of ellipticity with radius by 
stacking several galaxy samples. For cluster haloes, we can 
estimate the ellipticity as a function of radius for an indi- 
vidual galaxy cluster. We also point out that the minimum 
(^obs) provides another approach for determining the centre 
of galaxy clusters. Moreover, we introduced a weight func- 
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Figure 3. Tlic surface density k map of the simulated cluster 
taken from the PHOENIX simulations. This halo is used for gen- 
erating our flexion data. 
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annuli for the elliptical NFW halo (pluses) and the PHOENIX 
halo in 3 different projections, in units of arcsec^^. 



tion (in Appendix) to suppress the bias due to the asym- 
metric distribution of background galaxies. The variance of 
our estimate becomes signiflcantly smaller after employing 
the weight function. We derived the second order statistics 
of the flexion ratio (r^) and its analytical relation to (r). By 
combining (r) and {r^), in principle the bias due to intrinsic 
noise can be subtracted. We provided a method to reduce 
the bias based on a linear approximation (For more details, 
see Eq. AlO). 

We however notice that there are some difficulties in 
estimating the halo ellipticity. 

• First of all, the noise model that we employ in this paper 
includes only intrinsic noise and numerical noise in the simu- 
lations. In reality there are other sources of noise, e.g. noises 
due to the point spread function and pixelisation. Moreover, 
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Figure 4. Same as Fig. 1 for ttie mock galaxy cluster from the PHOENIX simulations projected along the x-axis (left),y-axis (middle) 
and z-axis (right). The diamonds (crosses) are the mean result over 50 realisations using halo without (with) subhaloes. 
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Figure 6. {r^^^) distribution due to the ccntroid offset for 10 
realisations. For each realisation, 100 random offsets are used. 



the noise of flexion measurement appears different using dif- 
ferent methods (Goldberg & Leonard 2007). Thus the noise 
in practice may be non-Gaussian. Therefore, the results pre- 
sented here are based on somewhat idealised observations. 

• The difficulty in reducing bias for small ellipticity haloes 
is mainly due to the lack of knowledge for the flexion noise 
behaviours. Moreover, the linear approximation that we em- 
ployed is based on the assumption that the bias is indepen- 
dent of the halo ellipticity and is smaller than the flexion 
ratio signal. In reality, neither of these assumptions may be 
valid. In particular, for small elliptical haloes, our measure- 
ment will be dominated by intrinsic and other noises. 

• The optimal binning to the data is also difficult since we 
do not know how the ellipticity changes with radius. A coarse 
binning can only estimate the "mean" ellipticity. Moreover, 
as we show in the figures, flexion data on different radii 
have different signal-to- noise ratio (S/N), thus better data 
at small radii may be polluted by the low S/N data at large 
radii. On the other hand, if many bins are used, not enough 
galaxies will be contained within each bin. 

• The degeneracy of substructures and halo ellipticity: 
substructures behave similarly to intrinsic flexion in our esti- 
mate, but at different magnitudes. Substructures are there- 
fore another kind of "noise". On the other hand, flexion 



can be used to measure substructures in galaxies or clusters 
(e.g. their mass function and spatial distribution). Fig. 5 is 
a rough estimate by comparing the ap magnitude without 
any knowledge of substructures or ellipticity. We can see 
that the combination of different approaches may constrain 
both shapes and substructures in dark matter haloes. Again 
good knowledge of flexion noise will be necessary. Also, a 
realistic model of the substructure mass function, and the 
spatial distribution of substructure will be helpful. 

A high number density of background galaxies is de- 
sirable. Current surveys such as COSMOS may have suffi- 
ciently high galaxy number densities to enable our method 
to constrain mass ellipticities on a few radial bins for an 
individual galaxy cluster. This is again based on a simple 
noise model. In reality, the noise of high magnitude images 
is larger. The number density of high quality images using 
current observations may not as high as we used. For galaxy- 
sized haloes, we can perform stacking gala:xy-galaxy lensing. 
Notice that different from the shear method, we do not need 
to align the major axis of the foreground galaxies. We do, 
however, assume that we can select dark matter haloes with 
similar shapes by carefully choosing the foreground galaxies 
according to their luminosity, colour and shape. The differ- 
ent morphologies will cause uncertainty in our estimate. On 
the other hand, the number density of galaxy-galaxy pairs 
will be very large in current and future surveys, i.e. the num- 
ber of stacked galaxies in reality can be much higher than 
what we used in our test. Fig. 2 indicates we can measure t 
to standard deviation Ot ~ 0.05 by stacking 20 foreground 
galaxies for e « 0.2. With a much larger number of galaxies 
available in future observations, the accuracy of this mea- 
surement can be improved. The improvement will be Pois- 
sonian if there are no other sources of systematic errors. An 
accurate measurement of how the ellipticity varies with ra- 
dius will provide a strong test of galaxy formation models. 
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APPENDIX A: 
REDUCTION 



BIAS ESTIMATE AND 



In this appendix, we will introduce two approaches that will 
reduce the estimated error. 



Al Weight function 

One of the biases comes from the asymmetric distribution 
of background galaxies. To account for such effects, we in- 
troduce a weight function 



1 



(Al) 



where A is a normalisation, and A''^ is the number of back- 
ground galax;ies in a polar angle range Aifi centred at the 
ith data. Therefore, the data is down-weighted where the 
number density is high. One can see that if the galaxy im- 
ages are circularly symmetrically distributed with respect 
to the foreground lens galaxy/cluster, the weight function 
will be uniform Wi = 1/A'^, where A'^ is the total number of 
background galaxies. The weight function is employed when 
calculating the expectation flexion ratio from data 



1=1 

iV 



(A2) 
(A3) 



A2 Bias reduction 

Due to the difficulty of measurement and our incomplete un- 
derstanding of the point spread function (PSF), our knowl- 
edge about the intrinsic flexion is not sufficient to give an 
accurate noise model. Since flexion has a dimension of in- 
verse length, the intrinsic flexion depends on the image size, 
and therefore the depth of survey. Moreover, Goldberg & 
Leonard (2007) found that different methods of measure- 
ment may lead to different noise behaviors. 

In ESll, it has been found that estimates using flexion 
ratios will be larger than the true ones due to the intrinsic 
noise. Here we adopt a simple model for the bias in the 
flexion ratio r, by assuming that the bias is independent of 
the halo ellipticity and lens strength. The observed mean 
flexion ratio can be written as 



(j-obs) = (r) + br, 



(A4) 



where br stands for the bias in the (r) estimate. Notice that 
br is not zero because the noise of the flexion ratio is not 
Gaussian (even if we assume that the noise of two flex- 
ion components is Gaussian, the ratio of two Gaussians no 
longer follows a Gaussian distribution), (r^) can be obtained 
through two ways. First, we can directly calculate (ro^s) 
from data using Eq.(A3). Secondly, according to Eqs.(6) and 
(8), we have 



(j-cal) 



(e 



T<r)/2 



2e'^('->/2 



(A5) 



We can thus calculate (r^^i) using Eq.(A5) from (robs)- 
The bias strength will be different in these two methods 
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Figure Al. Comparison of the ellipticity estimators 
with/without the weight function. In the top panel, the 
solid line is the input ellipticity, the pluses are the results using 
the weight function (Eq.Al), and the circles are the results 
without the weight function. The error bar is calculated by the 
standard deviation. The bottom panel represents the relative 
errors. 



(Fig. A2). One can see that both estimates give larger val- 
ues than the theoretical prediction (solid line), and they are 
different from each other. We make a flrst order expansion 

(rLs) = {r'} + 2{r)br + 0{bl); 



(rL) = {r'} + 



d{r) 



br + 0{bl 



ir) + -^W 



T{r)/2 



(A6) 
(A7) 

(A8) 



The difference between {rli,^) and (r^ai) provides an estimate 
of the bias br. 



A = (r-obs) - (r-cal) 

- (2(.) - ^ (c-«/^ 



(A9) 
-"<'->''2))6, + 0(&?). (AlO) 



We can subtract br from the estimate, and obtain a relatively 
clean result. However, this only works for large e. For small 
e (< 0.1), our approximation is no longer valid, since our 
measurements are dominated by noise. In the case of small 
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Figure A2. Comparison of {r"^^^) and {r^^i) with e: The solid 
line is for theoretical value of {r"^) (Eq.8), the pluses are obtained 
directly from data (Eq.A3), and the circles are calculated using 
Eq.(A5), with (r) obtained from data (Eq.A2). 

e, the h^ term is larger than the other terms in Eq.(A8), and 
thus cannot be neglected. Due to incomplete knowledge of 
the flexion noise, it is difficult to estimate when we should 
take the 6^ term into consideration. In practice, we first cal- 
culate (robs) and (robs) from the data, and then calculate 
{'^cai) using Eq.(A5). The bias is estimated through 



hr — 



Yobs) ~ Veal) 



2(robs) - f (e-<'-°b=>/2 - e— <'-°b=)/2) ' 



(All) 



If a negative br is obtained, we use for = "iKrobs) ^ (''cai)! as 
an empirical approximation. 



APPENDIX B: NUMERICAL TEST WITH A 
NON-SINGULAR ISOTHERMAL ELLIPSOID 
MODEL 

In this section we describe some simulations that we have 
performed in order to test the behaviours of the estimators 
given in the previous section. We model the halo surface 
mass density by a non-singular isothermal ellipsoid (NIE) 
profile: 



<e) 



(Bl) 



where the Einstein radius is 9e ~ 6 arcsec and the core 
radius is taken to he 8c = 2 arcsec. The redshifts of mock 
sources are infinite. We have 50 flexion data points in each 
realisation, which are randomly distributed on 1 x 1 arcmin^ 
"source plane". Data inside 6 arcsec and outside 30 arcsec 
are discarded (see Sect. 3.1). 

Bl Test of the weight function 

We generate 20 sets of mock data with different elliptici- 
ties e = 0.03i, i — 1,2, ...20. For each eUipticity, we use 50 
realisations. For each realisation, we calculate the expected 
flexion ratio in two ways. First, we calculate it directly using 
the mean of the data; second, we calculate (r) through our 



weight function (Eq.A2). Then we estimate e using Eq.(9). In 
Fig. Al, we show the mean estimate e over 50 realisations vs. 
the input values. The solid line shows the equality line, the 
pluses and circles are estimates with and without the weight 
function respectively. The error bars show the standard de- 
viations. Both estimates closely trace the input values, but 
the ones with weight function have smaller variances, es- 
pecially for large ellipticities. The bottom panel shows the 
relative error |5e/e|. One can see that the biases with the 
weight function are significantly smaller, and all are below 
2%. Hence throughout this paper the weight function has 
been employed. 



B2 The bias reduction 

We apply a simple model to generate noise in the data, 
jrobs = j-^ 4. ^^.j _|_ [(^jr^ _|_ jj^^) (as in ESll). The intrinsic 
flexion noises n/i,n/2 are drawn from a Gaussian distribu- 
tion independently, each component is characterized by a 
standard deviation aj^i = ctjf2 = 0.03 arcsec~^. 

We first compare (r^) using two methods. In Fig. A2, 
we show (r^bs) (pluses) and (r^ai) (circles) of different input 
e. For each e, we plot 5 realisations. The solid line is the 
theoretical prediction. Although the noise is generated in- 
dependent of e, the error and difference between (r^bs) and 
(r^jji) do increase with e. 

Next we perform a test on our approach to bias reduc- 
tion. 50 noise realisations for each e are generated. In each 
realisation, (robs) and (robs) are calculated directly from 
the mock data, and (r^ai) are calculated using Eq. (A5). 
In Fig. Bl we compare the estimated values with the input 
ones (solid line) . The pluses are the average of 50 realisations 
using Eq.(9). We also estimate the bias br using Eq.(AlO), 
and obtain an estimate from (r^ia^n) using Eq.(lO). We need 
to point out that for small e (in case of our noise model, 
e < 0.1), we adopt a modified approach since the linear 
approximation (Eq.AlO) is no longer valid. The results are 
shown by the circles in Fig. Bl. Overall, the bias is reduced, 
but for small e we still have larger errors. The bias reduction 
method is employed if not mentioned otherwise in Sect. 3. 
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